#include <iostream>
#include <vector>
#include <numeric>
#include <cmath>

using namespace std;

static vector<int> square_root_period(const int n) {
    vector<int> v;
    const int square_root = sqrt(n);
    int a = 1;
    int b = 0;
    int c = 1;
    while (true) {
        int r = (a * square_root + b) / c;
        v.push_back(r);
        if (r == 2 * square_root)
            return v;

        b -= c * r;
        int m = gcd(gcd(a, c), gcd(b, c));
        a /= m;
        b /= m;
        c /= m;

        int x = a * c;
        int y = -b * c;
        int z = a * a * n - b * b;
        a = x;
        b = y;
        c = z;
    }
    return {};
}

static pair<double, double> pell_equation(const int n) {
    vector<int> v = square_root_period(n);
    int k = v.size() - 1;
    for (int i = 1; i <= k; ++i)
        v.push_back(v[i]);

    double p0 = v[0];
    double q0 = 1;
    double p1 = v[1] * v[0] + 1;
    double q1 = v[1];
    k = k % 2 == 0 ? k : 2 * k;
    for (auto i = 2; i < k; ++i) {
        double p2 = v[i] * p1 + p0;
        double q2 = v[i] * q1 + q0;
        p0 = p1;
        q0 = q1;
        p1 = p2;
        q1 = q2;
    }

    return {p1, q1};
}

int main() {
    double max = 0; // if uint64_t, overflow 
    cout << "x^2 - D * y^2 = 1" << endl;
    for (int d = 2; d <= 1000; ++d) {
        if (sqrt(d) == floor(sqrt(d)))
            continue;

        auto [p, q] = pell_equation(d);
        if (p > max) {
            max = p;
            cout << "D = " << d << ", x = " << p << ", y = " << q << endl;
        }
    }
}
